class: title-slide, middle, right <br> <br> <br> .line_space_15[ ## .text_80[Modelos estructurales marginales para el control de<br>sesgos en estudios observacionales con factores de<br>riesgo y exposición tiempo-dependientes] ] <br> .line_space_11[ <br> .text_70[[Código en: <svg aria-hidden="true" role="img" viewBox="0 0 496 512" style="height:1em;width:0.97em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"/></svg>](https://github.com/AGSCL/DSPUCH)] .text_110[Seminario Métodos de Investigación en Salud Pública] ] .bg-text[ 10 de julio, 2022 .text_100[Andrés González Santa Cruz] .text_50[gonzalez.santacruz.andres@gmail.com] [<svg aria-hidden="true" role="img" viewBox="0 0 496 512" style="height:1em;width:0.97em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"/></svg>](https://github.com/AGSCL) [<svg aria-hidden="true" role="img" viewBox="0 0 512 512" style="height:1em;width:1em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:green;overflow:visible;position:relative;"><path d="M294.75 188.19h-45.92V342h47.47c67.62 0 83.12-51.34 83.12-76.91 0-41.64-26.54-76.9-84.67-76.9zM256 8C119 8 8 119 8 256s111 248 248 248 248-111 248-248S393 8 256 8zm-80.79 360.76h-29.84v-207.5h29.84zm-14.92-231.14a19.57 19.57 0 1 1 19.57-19.57 19.64 19.64 0 0 1-19.57 19.57zM300 369h-81V161.26h80.6c76.73 0 110.44 54.83 110.44 103.85C410 318.39 368.38 369 300 369z"/></svg>](https://orcid.org/0000-0002-5166-9121) ] <img src="data:image/png;base64,#./_style/Logo_nDP_monotono_vertical_en.png" width="15%" /> ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - Mi nombre es andrés gonzález y en esta instancia presentaré sobre los modelos estructurales marginales - Decidí abordar un tema metodológico no tratado en clases, que les será útil para mi tesis y que podría ser útil para otros. - --- layout: true class: animated, fadeIn --- ## Índice <br> 1. Problemática 2. Causalidad 3. Ajuste confusión 4. Marginal Structural Models 5. Aplicación 6. Desafíos ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA INDICE** *#_#_#_#_#_#_#_#_#_#_ Primeras dos, repaso: 1- **Origen de mi interés y problemática**: factores de riesgo y exposición tiempo-dependientes: Ejemplo de arbol de decisiones y que al final obtenía puros efectos condicionales. 2- **Causalidad**: La tarea de la epidemiología, definicion de causalidad y de confusión. Distribución marginal vs. conjunta, etc. Datos observacionales. Supuestos: SUTVA; Sesgos.Supuestos a la base, amenazas, sesgos, y particularmente sesgos que ocurren para este diseño de datos. 3- **Ajuste confusión**: Las formas. En una de esas, mostrar gif y animaciones (https://youtu.be/j8J2L_g76c4?t=68 o https://twitter.com/nickchk/status/1068215492458905600 / https://github.com/NickCH-K/causalgraphs/blob/master/Animation%20of%20IV.R / https://nickchk.com/causalgraphs.html). Explicar las diferencias entre estandarización y restricción --> ir poniendo el código paso por paso 4- **MSM**, los tipos.Explicar por qué son **estructurales** (1) y **marginales**(2), explicar las situaciones en que se ocupan y son útiles. Incluso con datos aleatorizados. Sesgos que resuelven 5- **Aplicación**: Explicar estructura hipotética de los datos, aplicación de un modelo simple. 6- **Desafíos y aproximaciones más complejas** --- ## Problemática .panelset.sideways[ .panel[.panel-name[Estructura hipotética más cercana a lo real] .details-code[ ```r # Libraries library(ggplot2) set.seed(2125) # Create data data <- data.frame( y=abs(rpois(1:250,15)), y2=abs(rpois(1:250,15)) ) %>% dplyr::filter(y2>y, y>=8, y2<=22) %>% dplyr::mutate(Paciente=row_number()) %>% #filtrar tratamientos más largos que 3 años dplyr::mutate(diff_treat=y2-y) %>% dplyr::filter(diff_treat<=3) for (i in 1:nrow(data)){ data$y3[i]<-base::sample(x=seq(from=data$y2[i]+1,to=22),1) data$y3[i]<-ifelse(data$y3[i]<=data$y2[i],22,data$y3[i]) data$y3[i]<-ifelse(data$y3[i]>=23,22,data$y3[i]) data$y4[i]<-ifelse(!is.na(data$y3[i]),base::sample(x=seq(from=data$y3[i]+1,to=22),1),22) data$y4[i]<-ifelse(data$y4[i]<=data$y3[i],22,data$y4[i]) data$y4[i]<-ifelse(data$y4[i]>=23,22,data$y4[i]) } set.seed(2125) pac_aleatorio1<-sample(1:max(data$Paciente),40) pac_aleatorio2<-setdiff(sample(1:max(data$Paciente),40), pac_aleatorio1) set.seed(2125) pac_aleatorio3<-unique(data$Paciente)[sample(1:length(data$Paciente),15)] # Horizontal version, antes era 1985 end_plot<-20 fig_trans<-ggplot(data) ###>=10 #datos de tratamiento en periodo de seguimiento, posterior al 2010 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y>=10,y2<=end_plot, Paciente %in% pac_aleatorio3), aes(x=Paciente, xend=Paciente, y=y, yend=y2), color="#EF9D2F", alpha=.6, size=.8)) fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y>=10,y2<=end_plot, !Paciente %in% pac_aleatorio3), aes(x=Paciente, xend=Paciente, y=y, yend=y2), color="#21177A", alpha=.6, size=.8)) #datos de tratamiento en periodo de seguimiento, exceden el seguiemiento. linea entera hasta el 2020 fig_trans<-try(fig_trans+ geom_segment(data=dplyr::filter(data,y>=10,y2>end_plot), aes(x=Paciente, xend=Paciente, y=y, yend=end_plot), color="#21177A", alpha=.6, size=.8)) #datos de tratamiento en periodo de seguimiento, exceden el seguiemiento. linea entrecortada despues del 2020 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y>=10,y2>end_plot), aes(x=Paciente, xend=Paciente, y=end_plot, yend=y2), color="#21177A", alpha=.6, linetype="dotted", size=.8)) ###<10 #datos de tratamiento en periodo de seguimiento, exceden el seguiemiento. linea entrecortada despues del 2020 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y<10), aes(x=Paciente, xend=Paciente, y=y, yend=10), color="#21177A", alpha=.6, linetype="dotted", size=.8)) #datos de tratamiento en periodo de seguimiento, menos de 2010 pero sólo y1 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y<10,y2>=10,y2<=end_plot, Paciente %in% pac_aleatorio3), aes(x=Paciente, xend=Paciente, y=10, yend=y2), color="#EF9D2F", alpha=.6, size=.8)) #fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y<10,y2>=10,y2<=end_plot, !Paciente %in% pac_aleatorio3), aes(x=Paciente, xend=Paciente, y=10, yend=y2), color="#21177A", alpha=.6, size=.8)) #datos de tratamiento en periodo de seguimiento, exceden el seguiemiento. linea entera hasta el 2020 # fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y<10,y2>end_plot), aes(x=Paciente, xend=Paciente, y=10, yend=end_plot), color="#21177A", alpha=.6, size=.8)) # CONTACTO: VICTIMARIO: debo saber si y3 está dentro o no del seguimiento ###>=10 #datos de justicia en periodo de seguimiento, linea solida fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y3<=end_plot, y4<=end_plot, Paciente%in% pac_aleatorio1), aes(x=Paciente, xend=Paciente, y=y3, yend=y4), color="#21177A", alpha=.6, size=.8)) #datos de justicia en periodo de seguimiento, exceden el seguiemiento pero parten en él. linea solida hasta el 2020 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y3<=end_plot, y4>end_plot, Paciente%in% pac_aleatorio1), aes(x=Paciente, xend=Paciente, y=y3, yend=end_plot), color="#21177A", alpha=.6, size=.8)) #si estoy dentro de los datos en y, peroen end plot no, tengo q hacer el interlineado afuera fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y3<=end_plot, y4>end_plot, Paciente%in% pac_aleatorio1), aes(x=Paciente, xend=Paciente, y=end_plot, yend=y4), color="#21177A", alpha=.6,linetype="dotted", size=.8)) #datos de tratamiento en periodo de seguimiento, exceden el seguiemiento. linea entrecortada despues del 2020 fig_trans<-try(fig_trans+geom_segment(data=dplyr::filter(data,y3>end_plot, y4>end_plot, Paciente%in% pac_aleatorio1), aes(x=Paciente, xend=Paciente, y=y3, yend=y4), color="#21177A", alpha=.6, linetype="dotted", size=.8)) fig_trans_final<- #Puntos fig_trans+ #Tratamientos completados geom_point(data=dplyr::filter(data,Paciente%in% sample(c(pac_aleatorio1,pac_aleatorio2),15), Paciente %in% pac_aleatorio3), aes(x=Paciente, y=y2), color="#EF9D2F", size=3, alpha=.6) + #geom_point(data=dplyr::filter(data,Paciente%in% pac_aleatorio1, Paciente %in% pac_aleatorio3), aes(x=Paciente, y=y4),color="#EF9D2F",size=3,alpha=.6)+ geom_point(data=dplyr::filter(data,Paciente%in% sample(c(pac_aleatorio1,pac_aleatorio2),15), !Paciente %in% pac_aleatorio3), aes(x=Paciente, y=y2), color="#21177A", size=3, alpha=.6) + geom_point(data=dplyr::filter(data,Paciente%in% pac_aleatorio1, !Paciente %in% pac_aleatorio3), aes(x=Paciente, y=y4),color="#21177A",size=3,alpha=.6)+ theme_light() + coord_flip() + #fechas donde yo tomo gente annotate("rect", xmin=-Inf, xmax=Inf, ymin=10, ymax=10.5, alpha = .2, fill="blue")+ annotate("rect", xmin=-Inf, xmax=Inf, ymin=19.5, ymax=20, alpha = .2, fill="blue")+ theme( panel.grid.major.y = element_blank(), panel.border = element_blank(), axis.ticks.y = element_blank() )+ scale_x_continuous(breaks=seq(1,max(data$Paciente),by=10))+ scale_y_continuous(breaks=seq(min(data$y),max(data$y4),by=2), labels=seq(min(data$y),max(data$y4),by=2)+2000)+ labs(y="Tiempo de seguimiento (en trimestres)", x="Individuos (ID)")+ #labs(y="Follow-up (in quarters)",x="Individual(ID)", caption="Note. Dot= Complete treatment;\nSquare= Contact w/ justice system (imputed);\nTriangle= Contact w/ justice system (victim);\nBlue line= Time in treatment;\nOrange line= Time in conctact w/justice system;\nShared area=Follow-up window")+ theme(plot.caption = element_text(hjust = 0, face= "italic"))+ labs(caption="Punto= Completa treatmento;\nAzul= Tratamiento ambulatorio;\nAmarillo= Tratamiento residencial; \nLínea= Tiempo en tratamiento;\nArea sombreada=Inicio y término seguimiento")+ theme( panel.background = element_rect(fill = "transparent", colour = NA_character_), # necessary to avoid drawing panel outline #panel.grid.major = element_blank(), # get rid of major grid #panel.grid.minor = element_blank(), # get rid of minor grid plot.background = element_rect(fill = "transparent", colour = NA_character_), # necessary to avoid drawing plot outline legend.background = element_rect(fill = "transparent"), legend.box.background = element_rect(fill = "transparent"), legend.key = element_rect(fill = "transparent")) invisible("Paciente 49: 2011 al 2014 cursa un primer tratamiento que abandona, luego del 2016 al 2018 completó un tratamiento") fig_trans_final ``` <img src="data:image/png;base64,#seminario_files/figure-html/p1-est-1.png" width="70%" height="70%" style="display: block; margin: auto;" /> ] ] - Pacientes cursando tratamiento por consumo problemático sustancias (2010-2019) - Usuarios pueden haber experimentado el evento más de una vez, en distintas oportunidades - Ej. Paciente N°55 completa tto. residencial y experimenta el evento de interés (readmisión) - Completar el tratamiento podría modificar .panel[.panel-name[Primera aproximación] Objetivo: **Estimar el efecto de la modalidad de tratamiento por trastornos por uso de sustancias (ambulatorios vs. residenciales) a la base en la probabilidad de experimentar subsiguientes readmisiones a TUS** Análisis propuestos en esa oportunidad: .superbigimage[ .text_50[ `$$readmisión\sim A_{1 (Modalidad\,base)}+X_{2 (Completa\,tto.\,base)}+\epsilon$$` `$$2^{da}\,readmisión\sim A_{1(Modalidad\,base)}+L_{2 (Completa\,tto.\,base)}+L_{3 (Completa\,2^{do}\,tto.)}+L_{4 (Días\,desde\,1^{er}\,ingreso)}+ \epsilon$$` `$$3^{era}\,readmisión\sim A_{1(Modalidad\,base)}+ L_{2 (Completa\,tto.\,base)}+ L_{3 (Completa\,2^{do}\,tto.)}+ L_{4 (Completa\,3^{er}\,tto.)}+ L_{5 (Días\,desde\,1^{er}\,ingreso)}+\epsilon$$` `$$4^{ta}\,readmisión\sim A_{1 (Modalidad\,base)}+ L_{2 (Completa\,tto.\,base)}+ L_{3 (Completa\,2^{do}\,tto.)}+ L_{4 (Completa\,3^{er}\,tto.)}+ L_{5 (Completa\,4^{to}\,tto.)}+ L_{6 (Días\,desde\,1^{er}\,ingreso)}+\epsilon$$` ] ] .center[**¿Cómo se reflejaría en un diagrama?**] ] .panel[.panel-name[Árbol causal] <img src="data:image/png;base64,#seminario_files/figure-html/unnamed-chunk-2-1.png" width="550" style="display: block; margin: auto;" /> ] .panel[.panel-name[Esquema Diagrama Causal Simplificado] .details-code[ ```r dag23 <- dagitty('dag { bb="0,0,1,1" "A0 (Modalidad tto. basal)" [exposure,pos="0.184,0.608"] "A1 (Modalidad tto. 2)" [pos="0.451,0.605"] "A2 (Modalidad tto. 3)" [pos="0.734,0.605"] "L0 (Vector car. individuales)- Pareamiento" [adjusted,pos="0.121,0.397"] "L1 (Días previos a la readmisión 2)" [adjusted,pos="0.507,0.341"] "L2 (Días previos a la readmisión 3)" [adjusted,pos="0.691,0.365"] "LM0(Completa tto. basal)" [adjusted,pos="0.269,0.707"] "LM1 (completa tto. 2)" [adjusted,pos="0.566,0.714"] "LM2 (completa tto. 3)" [adjusted,pos="0.838,0.722"] "U0 (ej., factores genéticos)" [latent,pos="0.124,0.273"] "U1 Car. indivudales (ej., patrón de consumo)" [latent,pos="0.444,0.829"] "U2 Car. individuales T2" [latent,pos="0.733,0.814"] "Y1(Readmisión)" [outcome,pos="0.355,0.604"] "Y2(2da Readmisión)" [outcome,pos="0.626,0.605"] "Y3(4ta Readmisión)" [outcome,pos="0.960,0.607"] "A0 (Modalidad tto. basal)" -> "LM0(Completa tto. basal)" "A0 (Modalidad tto. basal)" -> "Y3(4ta Readmisión)" "A1 (Modalidad tto. 2)" -> "LM1 (completa tto. 2)" "A2 (Modalidad tto. 3)" -> "LM2 (completa tto. 3)" "L0 (Vector car. individuales)- Pareamiento" -> "A0 (Modalidad tto. basal)" "L0 (Vector car. individuales)- Pareamiento" -> "Y1(Readmisión)" "L1 (Días previos a la readmisión 2)" -> "Y2(2da Readmisión)" "L2 (Días previos a la readmisión 3)" -> "Y3(4ta Readmisión)" "LM0(Completa tto. basal)" -> "U1 Car. indivudales (ej., patrón de consumo)" "LM0(Completa tto. basal)" -> "Y1(Readmisión)" "LM0(Completa tto. basal)" -> "Y2(2da Readmisión)" "LM0(Completa tto. basal)" -> "Y3(4ta Readmisión)" "LM1 (completa tto. 2)" -> "U2 Car. individuales T2" "LM1 (completa tto. 2)" -> "Y2(2da Readmisión)" "LM1 (completa tto. 2)" -> "Y3(4ta Readmisión)" "LM2 (completa tto. 3)" -> "Y3(4ta Readmisión)" "U0 (ej., factores genéticos)" -> "L0 (Vector car. individuales)- Pareamiento" "U0 (ej., factores genéticos)" -> "Y1(Readmisión)" "U0 (ej., factores genéticos)" -> "Y2(2da Readmisión)" "U0 (ej., factores genéticos)" -> "Y3(4ta Readmisión)" "U1 Car. indivudales (ej., patrón de consumo)" -> "A1 (Modalidad tto. 2)" "U1 Car. indivudales (ej., patrón de consumo)" -> "LM1 (completa tto. 2)" "U2 Car. individuales T2" -> "A2 (Modalidad tto. 3)" "U2 Car. individuales T2" -> "LM2 (completa tto. 3)" "Y1(Readmisión)" -> "A1 (Modalidad tto. 2)" "Y1(Readmisión)" -> "U1 Car. indivudales (ej., patrón de consumo)" "Y2(2da Readmisión)" -> "A2 (Modalidad tto. 3)" "Y2(2da Readmisión)" -> "U2 Car. individuales T2" }') tidy_dag23 <- tidy_dagitty(dag23) %>% dplyr::mutate(label=dplyr::case_when(grepl("A0",as.character(name))~"A0", grepl("A1",as.character(name))~"A1", grepl("A2",as.character(name))~"A2", grepl("L0",as.character(name))~"L0", grepl("L1",as.character(name))~"L1", grepl("L2",as.character(name))~"L2", grepl("LM0",as.character(name))~"LM0", grepl("LM1",as.character(name))~"LM1", grepl("LM2",as.character(name))~"LM2", grepl("U0",as.character(name))~"U0", grepl("U1",as.character(name))~"U1", grepl("U2",as.character(name))~"U2", grepl("Y1",as.character(name))~"Y1", grepl("Y2",as.character(name))~"Y2", grepl("Y3",as.character(name))~"Y3", T~as.character(name))) %>% dplyr::mutate(label2=dplyr::case_when(grepl("U",name)~"latent",grepl("L|LM",name)~"adj",grepl("LM",name)~"white", T~"black")) %>% dplyr::mutate(adjusted=factor(dplyr::case_when(grepl("L|LM",name)~"adjusted",T~"unadjusted"))) edge_function <- ggdag:::edge_type_switch("link_arc") dag23_plot<- ggdag:::if_not_tidy_daggity(tidy_dag23) %>% ggdag:::node_status() %>% ggplot2::ggplot(ggplot2::aes(x = x, y = y, xend = xend, yend = yend, color = status, shape=factor(adjusted)))+ #edge_function()+ scale_adjusted()+ ggdag:::breaks(c("exposure", "outcome","latent"))+ geom_dag_edges_arc(curvature = c(0,5,rep(0,26)))+ ggdag:::geom_dag_point(size = 16)+ ggdag:::geom_dag_label_repel(ggplot2::aes_string(label = "label", fill = "status"), size = 4.88, col = "white", show.legend = FALSE)+ theme_dag()+ scale_shape_manual(values = c(15, 16), name="Ajustado", labels=c("Sí", "No"))+ scale_fill_manual(values = c("#003891", "#EF9D2F","gray30"), name="Estatus",na.value="black", labels=c("Exposición", "Resultado","No observado"), limits = c('exposure', 'outcome','latent'))+ scale_color_manual(values = c("#003891", "#EF9D2F","gray30"), name="Estatus",na.value="black", labels=c("Exposición", "Resultado","No observado"), limits = c('exposure', 'outcome','latent'))+#E6E6E6 guides(linetype="none", edge_alpha="none", shape="none")+ guides(color=guide_legend(override.aes = list(arrow = NULL)))+#,guide_colourbar(order = 1) theme(plot.caption = element_text(hjust = 0))+ theme(legend.position = "bottom", aspect.ratio=6/10)+ labs(caption="Nota. Ak= Modalidad (Residencial/Ambulatoria); LMk= Completa tratamiento;\nL0= Confusores, características individuales y del centro a la base;Lk+= Días previos a la readmisión;\nYk= Resultadotratamiento; Uk=Car. no observadas para cada tiempo") dag23_plot ``` <img src="data:image/png;base64,#seminario_files/figure-html/p2-dag-1.png" height="450" style="display: block; margin: auto;" /> ] <div class="center"> ¿Qué ocurre con las modalidades en los tratamientos sucesivos?, ¿Confunden la asociación entre `\(A_0\)` ` e` `\(Y_3\)` `? `, ¿Sigue siendo causal en presencia de { `\(A_1\,, A_2\,,A_3\)` `}? </div> <!--- U puede tener relación con derivaciones no observadas posterior al tratamiento, que sería como un U1.5---> <!--- U puede tener relación con derivaciones no observadas posterior al tratamiento, que sería como un U1.5---> ] .panel[.panel-name[Entonces...] - Varias puertas traseras sin cerrar, confusión residual, abandono/pérdida, sobreajuste, etc. - El diseño y la estrategia analítica hasta el momento no permite responder a una pregunta causal de manera adecuada - Debiesen analizarse otras estrategias que capturen efectos longitudinales <br> .pull_r_30[ .text_60[ ***"In longitudinal studies with time-dependent confounding, identifying the structure allows us to detect situations in which stratification-based methods would adjust for confounding at the expense of introducing selection bias (p. 622)"*** `\(^{[1]}\)` ] ] <br> <br> <br> .down_center[ <div class="red"> ¿Se puede responder a la pregunta de si haber asistido a tratamiento residencial a la base ( `\(A_0\)` `) es beneficioso?, y si lo es, ¿qué régimen/estrategia es óptima o casi óptima? `\(^{[2]}\)` </div> ] ] ] ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - Causal tree graph, L=0 y L=1 son los resultados del tratamiento (completa/no-completa), y los A son residencial/ambulatorio. Y= es readmisión - La figura muestra distintas trayectorias de la cohorte. La entrada de paciente a la cohorte retrospectiva empieza a la fecha en que es admitido por primera vez en un tratamiento por uso de sustancias en el período de 2010-2019 y figura en las listas anuales del SENDA al 2010 al 2019 (independientemente si tienen tratamientos previos). - **Censura administrativa**: entrega información SENDA (2019-11-13), - Otra censura ocurre después que se produce un evento de resultado competidor - Cuando un paciente deja la cohorte sin experimentar resultados (ej, si se va fuera del país). - **Eventuales competidores**: tiempo en tratamiento (especialmente residenciales) o una sentencia **aunque debemos ver si es posible conocer este tiempo de censura**, muerte #_#_#_#_#_#_#_#_#_#_#_#_#_ The total effect cannot be estimated due to adjustment for an intermediate or a descendant of an intermediate. - El matching eElimina puertas traseras pero no se puede graficar. - **Hasta el momento es encesario hacer modificaciones al diseño y la estrategia analítica** --- class: center, middle ## Objetivo Este proyecto apunta a servir como un material introductorio sobre modelos estructurales marginales para el control de sesgos en estudios observacionales con factores de riesgo y/o exposición tiempo-dependientes. -- - Se mostrará su estimación paso a paso y mediante un ejercicio reproducible. -- - Se utilizará una base de datos simulada para ilustrar el proceso. ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ Ejemplos de publicaciones de este tipo hay: - Tutoriales de herramientas analíticas innovadoras en "Psychological Methods" (APA) https://www.apa.org/pubs/journals/met/call-for-papers-tutorials - Bradburn, M., Clark, T., Love, S. et al. Survival Analysis Part I: Multivariate data analysis – an introduction to concepts and methods. British Journal of Cancer 89, 431–436 (2003). https://doi.org/10.1038/sj.bjc.6601119 - Zhang, Z., Abarda, A., Contractor, A., Wang, J., & Dayton, C. (2018). Exploring heterogeneity in clinical trials with latent class analysis. Annals Of Translational Medicine, 6(7), 9. doi:10.21037/atm.2018.01.24 --- ## Causalidad (1) .pull_left[ - Salud pública- epidemiología `\(^{[3]}\)` - Qué pasaría si... // Qué causa ... `\(^{[4; 5]}\)` <img src="data:image/png;base64,#seminario_files/figure-html/unnamed-chunk-3-1.png" width="350" style="display: block; margin: auto auto auto 0;" /> - Asociación vs. Causalidad `\(^{[4; 6]}\)` - Inferencia `\(^{[6]}\)` - RCTs: gold-standard `\(^{[7]}\)` - Estudios observacionales .text_70[ ***"It’s much easier to get a result than it is to get an answer"*** `\(^{[8]}\)` ] ] .pull_right[ . <img src="data:image/png;base64,#seminario_files/figure-html/unnamed-chunk-4-1.png" width="350" style="display: block; margin: auto 0 auto auto;" /> .line_space_05[ .text_60[ .footnote[ Riseberg,E., Melamed, R., James, K., Alderete, T. & Corlin, L. (2021). Development and application of an evidence-based directed acyclic graph to evaluate the associations between metal mixtures and cardiometabolic outcomes. doi: https://doi.org/10.1101/2021.03.05.21252993 ] ] ] ] ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** (1) *#_#_#_#_#_#_#_#_#_#_ - Distribución marginal vs. conjunta - (a veces, no factibles, poco éticos, potencia estadística, etc.) - estandarizar el riesgo de la distribución conjunta (joint distribution) de los confusores en todo el seguimiento densidad conjunta de los datos observados para generar resultados potenciales bajo distintos escenarios de resultados --- ## Causalidad (2) **SUTVA** `\(^{[9; 10; 11]}\)` <img src="data:image/png;base64,#./_figs/sutva.jpg" width="300" style="display: block; margin: auto auto auto 0;" /> -- - *Consistencia* -- - *Intercambiabilidad* -- - *Positividad* ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** (2) *#_#_#_#_#_#_#_#_#_#_ - Si queremos que las estimaciones de un efecto puedan tener una interpretación causal, es necesario asumir algunos supuestos sobre los datos y su proceso de generación. potential outcomes framework (Neyman, 1923; Rubin, 1974). AT E = E[Y (1) -Y(0)] - estimates of effect that may fail to have a causal interpretation, even when confounding by unmeasured factors and model specification are absent - Asumimos NO Measurement Bias - Explicar causal approach, contrafactual, ajuste por confusión, hay métodos que restringen y otros que estandarizan - En las clases abordamos a structural approach to selection bias y en What if se mencionó someramente sobre estos mecanismos. - De ahí que surgieron algunas dudas --- ## Ajuste confusión - **Estratificación / Regresiones** - **Estandarización** `\(^{[12; 5; 13]}\)` ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ **NO CREO NECESARIO GRAFICAR LA DIFERENCIA ENTRE UNO Y OTRO (AL FINAL FUE MUY DIFÍCIL, PORQUE LOS EJEMPLOS NO ERAN SUFICIENTEMENTE BUENOS** - El ajuste de confusión es una forma de cerrar puertas traseras, siempre y cuando se cumplan otros supuestos (no hay error de medición, variables no-observadas (confusión residual), mediadores afectad) - Cuando se ajusta por variables en un modelo de regresión se obtiene un efecto condicional. El efecto manteniendo constante el resto de las variables. - The conditional approach handles confounders using stratification or modeling (e.g., adding covariates to be regressed to the outcome). --- ## Marginal structural models - ¿Por qué el nombre? - Tres principales: - IPTW - Fórmula G paramétrica - Estimación G ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - Que no haya Ajuste por variables intermedias o descendientes de una variable intermedia ("the mantra not to control for factors affected by exposure".- Overadjustment)Schisterman2009 - Control for `W` by seeing what `W` explains (perhaps with a regression) and taking it out - *crossover effects --- ## Aplicación (1) Se generó una base de datos de 1.000 observaciones con las siguientes características `\(^{[14]}\)` : - t1 = `Tratamiento residencial Tiempo 0 ` ( `\(A_0\)` ) ` ∼ Bernoulli(p=.23) - t2 = `Tratamiento residencial Tiempo 1 ` ( `\(A_1\)` ) ` ∼ Bernoulli(p=.22) - t3 = `Tratamiento residencial Tiempo 2 ` ( `\(A_2\)` ) ` ∼ Bernoulli(p=.15) - c1 = `Policonsumo` ∼ Bernoulli(p=.32) - c2 = `Edad de inicio consumo de sustancias ` ∼ N(3,0.5) - dt1 = `Meses sin readmisión Primer tto.` ∼ Poisson(λ=4) - dt2 = `Meses sin readmisión Segundo tto.` ∼ Poisson(λ=3) - dt3 = `Meses sin readmisión Tercer tto.` ∼ Poisson(λ=2) - v1 = `Confusor tiempo-dependiente (t1)` ∼ `\(5+ 0.4\beta_{Tto\,residencial\,basal} +0.78\beta_{Meses\,sin\,readmision\,1}+ N(0, \sqrt{0.99})\)` - v2 = `Confusor tiempo-dependiente (t2)` ∼ `\(5+ 0.4\beta_{Tto\,residencial\,(2^{do})} +0.78\beta_{Meses\,sin\,readmision\,1}+ N(0, \sqrt{0.55})\)` - v3 = `Confusor tiempo-dependiente (t3)` ∼ `\(5+ 0.4\beta_{Tto\,residencial\,(3^{er})} +0.78\beta_{Meses\,sin\,readmision\,1}+ N(0, \sqrt{0.33})\)` ` .details-code[ ```r #Semilla para replicar los resultados set.seed(2125) #definir confusor c1 (policonsumo al ingreso==1) policonsumo <-rbinom(1000,1,0.32) #definir confusor c2 (edad de inicio de consumo) edad_ini <-round(exp(rnorm(1000, 3, 0.5))) #definir tratamiento residencial al tiempo 0 t1 <-rbinom(1000,1,0.23) #definir tratamiento residencial al tiempo 1 t2 <-rbinom(1000,1,0.22) #definir tratamiento residencial al tiempo 2 t3 <-rbinom(1000,1,0.15) #meses libre de readmisión T1 dt1 <-rpois(1000,4) #meses libre de readmisión T2 dt2 <-rpois(1000,3) #meses libre de readmisión T3 dt3 <-rpois(1000,1) #definir el confusor tiempo-dependiente v1 como una función de t1 y dt1 v1 <- (0.4*t1 + 0.78*dt1 + rnorm(1000, 0, sqrt(0.99))) + 5 #definir el confusor tiempo-dependiente v2 como una función de t2 y dt2 v2 <- (0.4*t2 + 0.78*dt2 + rnorm(1e3, 0, sqrt(0.55))) + 5 #definir el confusor tiempo-dependiente v3 como una función de t3 y dt3 v3 <- (0.4*t3 + 0.78*dt3 + rnorm(1e3, 0, sqrt(0.33))) + 5 datos <- data.frame(policonsumo, edad_ini, v1, v2, v3, t1, t2, t3, dt1, dt2, dt3) %>% # Generamos una columna id para cada sujeto dplyr::mutate(id = row_number()) %>% # Convertimos la base de datos desde formato ancho a largo de todas las variables tiempo-dependientes. Aclaramos que la estructura de los nombres lleva inscrito primero el valor de interés, seguido por el tiempo/ola de aplicación. Y obtenemos el tiempo tomando el valor dígito de cada columna tidyr::pivot_longer(cols=c("v1", "v2", "v3", "t1", "t2", "t3", "dt1", "dt2", "dt3"), names_to = c('.value',"time"), names_pattern = "(.*)(\\d+)$") %>% # Ordenamos la base de datos por id y tiempo/ola dplyr::arrange(id, time) %>% # reordenamos las columnas dplyr::select(id, time, policonsumo, edad_ini, t, v, dt) %>% # formateamos la base de datos como data.frame as.data.frame() %>% dplyr::mutate(time=as.numeric(time)) #Primeros 4 filas head(datos,4) %>% knitr::kable("html", caption="Muestra de la base de datos (Primeras 4 filas)") %>% kableExtra::kable_classic(bootstrap_options = c("striped", "hover"),font_size = 10) ``` <table class=" lightable-classic" style='font-size: 10px; font-family: "Arial Narrow", "Source Sans Pro", sans-serif; margin-left: auto; margin-right: auto;'> <caption style="font-size: initial !important;">Muestra de la base de datos (Primeras 4 filas)</caption> <thead> <tr> <th style="text-align:right;"> id </th> <th style="text-align:right;"> time </th> <th style="text-align:right;"> policonsumo </th> <th style="text-align:right;"> edad_ini </th> <th style="text-align:right;"> t </th> <th style="text-align:right;"> v </th> <th style="text-align:right;"> dt </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 30 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 10.633698 </td> <td style="text-align:right;"> 7 </td> </tr> <tr> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 30 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 6.549493 </td> <td style="text-align:right;"> 3 </td> </tr> <tr> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 3 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 30 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 7.595415 </td> <td style="text-align:right;"> 3 </td> </tr> <tr> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 44 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 9.276188 </td> <td style="text-align:right;"> 2 </td> </tr> </tbody> </table> ] ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - Se asume que todos cursaron el tratamiento y no hay censura, que la velocidad es constante y los riesgos también. - Asumiendo que los riesgos son constantes y proporcionales (Pierce) - More generally, assuming constant hazard rates over fixed time intervals (known as a piecewise-exponential model) you can fit fairly flexible survival models in the form of poisson GLM --- ## Aplicación (2) .details-code[ ```r #```{r, echo=F, dev.args = list(bg = 'transparent'), fig.align="center", out.width=550, error=T} #fig.width = 7, fig.height = 5 muy grande dag34 <- dagitty('dag { bb="0,0,1,1" L1_policonsumo [pos="0.117,0.205"] L2_edad_ini [pos="0.222,0.112"] dt1 [outcome,pos="0.294,0.491"] dt2 [outcome,pos="0.572,0.491"] dt3 [outcome,pos="0.838,0.491"] t1 [exposure,pos="0.092,0.491"] t2 [exposure,pos="0.380,0.491"] t3 [exposure,pos="0.657,0.491"] v1 [pos="0.289,0.686"] v2 [pos="0.576,0.668"] v3 [pos="0.849,0.657"] L1_policonsumo -> dt1 L1_policonsumo -> dt2 L1_policonsumo -> dt3 L1_policonsumo -> t1 L1_policonsumo -> t2 L1_policonsumo -> t3 L2_edad_ini -> dt1 L2_edad_ini -> dt2 L2_edad_ini -> dt3 L2_edad_ini -> t1 L2_edad_ini -> t2 L2_edad_ini -> t3 t1 -> dt1 t1 -> t2 [pos="0.210,0.652"] t2 -> dt2 t2 -> t3 [pos="0.514,0.677"] t3 -> dt3 v1 -> dt1 v1 -> t1 v1 -> t2 v2 -> dt2 v2 -> t2 v2 -> t3 v3 -> dt3 v3 -> t3 }') tidy_dag34 <- tidy_dagitty(dag34) %>% dplyr::mutate(label=dplyr::case_when(grepl("A0",as.character(name))~"A0", T~as.character(name))) %>% dplyr::mutate(label2=dplyr::case_when(grepl("L|v",name)~"adj",grepl("LM",name)~"white", T~"black")) %>% dplyr::mutate(adjusted=factor(dplyr::case_when(grepl("L|v",name)~"adjusted",T~"unadjusted"))) edge_function <- ggdag:::edge_type_switch("link_arc") dag34_plot<- ggdag:::if_not_tidy_daggity(tidy_dag34) %>% ggdag:::node_status() %>% ggplot2::ggplot(ggplot2::aes(x = x, y = y, xend = xend, yend = yend, color = status, shape=factor(adjusted)))+ #edge_function()+ scale_adjusted()+ ggdag:::breaks(c("exposure", "outcome","latent"))+ theme_dag()+ geom_dag_edges_arc(curvature = c(rep(0,13),5,0,5,rep(0,9)))+ #14 y 16 de 24 ggdag:::geom_dag_point(size = 16)+ ggdag:::geom_dag_label_repel(ggplot2::aes_string(label = "label", fill = "status"), size = 4.88, col = "white", show.legend = FALSE)+ scale_shape_manual(values = c(15, 16), name="Ajustado", labels=c("Sí", "No"))+ scale_fill_manual(values = c("#003891", "#EF9D2F","gray30"), name="Estatus",na.value="black", labels=c("Exposición", "Resultado"), limits = c('exposure', 'outcome'))+ scale_color_manual(values = c("#003891", "#EF9D2F","gray30"), name="Estatus",na.value="black", labels=c("Exposición", "Resultado"), limits = c('exposure', 'outcome'))+#E6E6E6 guides(linetype="none", edge_alpha="none", shape="none")+ guides(color=guide_legend(override.aes = list(arrow = NULL)))+#,guide_colourbar(order = 1) theme(plot.caption = element_text(hjust = 0))+ theme(legend.position = "bottom", aspect.ratio=6/10)+ labs(caption="Nota. Ak= Modalidad (Residencial/Ambulatoria); dtk= Meses libre de readmisión;\nvk= Confusor tiempo-dependiente") dag34_plot ``` <img src="data:image/png;base64,#seminario_files/figure-html/tab-comp0_dag-1.png" height="450" style="display: block; margin: auto;" /> ] ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - E[Treatment=1, time=1 | Treatment=1, time=(time-1)] - Auto-regressive (AR1) correlation structure ... correlation to decay as the outcome values were farther away from the time of interest. --- ## Aplicación (3) . - Se estima una ponderación inversa al tratamiento que tiene en cuenta exposición y/o confusores tiempo-dependiente `\(^{[15]}\)`. .details-code[ ```r #Paquetes necesarios para generar los análisis library(geepack) library(survey) library(ipw) library(reshape) library(MuMIn) ########################################################################################### #Estimate ipw weights (time-varying) ########################################################################################### #Estimate inverse probability weights to fit marginal structural models, with a time-varying exposure and time-varying confounders. Within each unit under observation this function computes inverse probability weights at each time point during follow-up. The exposure can be binomial, multinomial, ordinal or continuous. Both stabilized and unstabilized weights can be estimated. # estimate inverse probability weights (time-varying) using a logistic regression #https://search.r-project.org/CRAN/refmans/ipw/html/ipwtm.html w <- ipwtm( exposure = t, family = "binomial", link = "logit", # Time invariant stuff numerator = ~ factor(policonsumo) + edad_ini, # All confounders denominator = ~ v + factor(policonsumo) + edad_ini, id = id, timevar=time, type="first", data = datos) summary(w$ipw.weights) #Se incorporan los pesos a la base de datos datos<- dplyr::bind_cols(datos, ipw=w$ipw.weights) ########################################################################################### # Para obtener las salidas de los coeficientes ########################################################################################### summary.geem <- function(object, ...) { Coefs <- matrix(NA,nrow=length(object$beta),ncol=5) Coefs[,1] <- c(object$beta) naive <- is.character(object$var) if(!naive && any(diag(object$var) < 0) ){ naive <- TRUE warning("Some elements of robust variance estimate < 0. Reporting model based SE.") } Coefs[,2] <- sqrt(diag(object$naiv.var)) if(naive){Coefs[,3] <- rep(NA, length(object$beta))}else{Coefs[,3] <- sqrt(diag(object$var))} if(naive){Coefs[,4] <- Coefs[,1]/Coefs[,2]}else{Coefs[,4] <- Coefs[,1]/Coefs[,3]} Coefs[,5] <- round(2*pnorm(abs(Coefs[,4]), lower.tail=F), digits=8) colnames(Coefs) <- c("Estimates","Model SE","Robust SE", "wald", "p") summ <<- list(beta = Coefs[,1], se.model = Coefs[,2], se.robust = Coefs[,3], wald.test = Coefs[,4], p = Coefs[,5], alpha = object$alpha, corr = object$corr, phi = object$phi, niter = object$niter, clusz = object$clusz, coefnames = object$coefnames, weights=object$weights, biggest.R.alpha = object$biggest.R.alpha) class(summ) <- 'summary.geem' return(summ) } ########################################################################################### # Primer modelo GEE - GEE con errores estándar robustos (no IPTW) ########################################################################################### mod1<-geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id ,data = data.frame(datos), family=poisson, corstr="ar1") summary(mod1) summ_mod1<- cbind.data.frame(cc=summ$coefnames, Estimate=summ$beta,Std.err=summ$se.robust) %>% dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>% dplyr::mutate(across(where(is.numeric),~round(exp(.),2))) invisible("Quasipoisson") mod2<-geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id, data = datos, family = MASS::negative.binomial(2), corstr = "ar1") #The known value of the additional parameter, theta. mod3<-geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id, data = datos, family = MASS::negative.binomial(1), corstr = "ar1") mod4<-geepack::geese(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id, data = datos, corstr = "ar1", family=poisson) ``` ] - Una vez que se genera el modelo de ecuaciones de estimación generalizadas (GEE) sin considerar la probabilidad inversa de asignación a tratamiento, se contrasta con el modelo en el que se pondera por dicha probabilidad inversa. Se seleccionó el modelo con menor QIC. .details-code[ ```r #https://gist.github.com/avallecam/56af06f46e5544c3af0f46344df20989 #https://stackoverflow.com/questions/49089476/any-updates-to-model-negative-binomial-distribution-data-with-gee-in-r #https://www.rdocumentation.org/packages/geepack/versions/1.3.4/topics/geese #https://stackoverflow.com/questions/13946540/negative-binomial-in-gee?rq=1 mod1b<- geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id ,weights=ipw, data = data.frame(datos), family=poisson, corstr="ar1") summary.geem(mod1b) summ_mod1b<- cbind.data.frame(cc=summ$coefnames, Estimate=summ$beta,Std.err=summ$se.robust) %>% dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>% dplyr::mutate(across(where(is.numeric),~round(exp(.),2))) mod2b<- geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id, weights=ipw, data = datos, family = MASS::negative.binomial(2), corstr = "ar1") #The known value of the additional parameter, theta. summary(mod2b) summ_mod2b<- cbind.data.frame(cc=summ$coefnames, Estimate=summ$beta,Std.err=summ$se.robust) %>% dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>% dplyr::mutate(across(where(is.numeric),~round(exp(.),2))) mod3b<- geem(dt~ t+ time+ factor(policonsumo)+ edad_ini, id=id, weights=ipw, data = datos, family = MASS::negative.binomial(1), corstr = "ar1") summary(mod3b) summ_mod3b<- cbind.data.frame(cc=summ$coefnames, Estimate=summ$beta,Std.err=summ$se.robust) %>% dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>% dplyr::mutate(across(where(is.numeric),~round(exp(.),2))) #You can use the quasi-likelihood under the independence model criterion (QIC; Pan, 2001) to assess the model fit. model.sel(mod1b, mod2b, mod3b, rank = QIC) ``` ] Por último, una tabla comparativa de las estimaciones del tratamiento. .details-code[ ```r mt<- rbind.data.frame( summ_mod1[which(summ_mod1$cc=="t"),], summ_mod1b[which(summ_mod1$cc=="t"),] ) rownames(mt) <- c("Modelo sin IPTW", "Modelo con IPTW") mt %>% knitr::kable("html", caption="Comparación resultados") %>% kableExtra::kable_classic(bootstrap_options = c("striped", "hover"),font_size = 10) ``` <table class=" lightable-classic" style='font-size: 10px; font-family: "Arial Narrow", "Source Sans Pro", sans-serif; margin-left: auto; margin-right: auto;'> <caption style="font-size: initial !important;">Comparación resultados</caption> <thead> <tr> <th style="text-align:left;"> </th> <th style="text-align:left;"> cc </th> <th style="text-align:right;"> Estimate </th> <th style="text-align:right;"> Std.err </th> <th style="text-align:right;"> lwr </th> <th style="text-align:right;"> upr </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> Modelo sin IPTW </td> <td style="text-align:left;"> t </td> <td style="text-align:right;"> 0.99 </td> <td style="text-align:right;"> 1.03 </td> <td style="text-align:right;"> 0.94 </td> <td style="text-align:right;"> 1.05 </td> </tr> <tr> <td style="text-align:left;"> Modelo con IPTW </td> <td style="text-align:left;"> t </td> <td style="text-align:right;"> 0.88 </td> <td style="text-align:right;"> 1.03 </td> <td style="text-align:right;"> 0.83 </td> <td style="text-align:right;"> 0.94 </td> </tr> </tbody> </table> ] ??? *#_#_#_#_#_#_#_#_#_#_ **NOTA** *#_#_#_#_#_#_#_#_#_#_ - applying inverse probability of treatment weights (IPTW) to the observations, which are then fitted to a marginal structural model (MSM). We simulate a simple dataset with a time-varying exposure. The data, collected at two time points, t = (0, 1), include a binary outcome measured at the end of follow up Y (1), a binary exposure at two time points {X(0),X(1)}, and a binary confounder of the exposure-outcome relationship that is also affected by prior exposure L(1). L(1) is related to the outcome by the common cause U, which we consider to be unmeasured. This particular simulation setup has been used previously to illustrate other causal inference methods.(Robins y Hernan, 2009) It addresses a central feature of the Bayesian g-formula: control of confounding by covariates that may also be affected by exposure. --- ## Desafíos - Structural Nested Models - Confusores tiempo-modificados (*time-modified confounding*) `\(^{[16]}\)` - Generar modelos estructurales de análisis de supervivencia (incluir probabilidad inversa de censura) - Variaciones - Modelos doble o triplemente robustos - Machine Learning (bagging, boosting, random forests, neural networks) --- class: center, middle # Gracias! <br> <div class="centered"> Contacto: gonzalez.santacruz.andres@gmail.com </div> <br> <br> <br> .down_right[ <img src="data:image/png;base64,#seminario_files/figure-html/unnamed-chunk-6-1.png" width="200" style="display: block; margin: auto;" /> ] .down_left[ <img src="data:image/png;base64,#./_style/Logo_nDP_color_hz_en.png" width="200" style="display: block; margin: auto;" /> ] --- ## Referencias . [1] M. A. Hernán, S. Hernández-Díaz, and J. M. Robins. "A structural approach to selection bias.". Eng. In: _Epidemiology (Cambridge, Mass.)_ 15.5 (sept.. 2004), pp. 615-625. ISSN: 1044-3983 (Print). DOI: [10.1097/01.ede.0000135174.63482.43](https://doi.org/10.1097%2F01.ede.0000135174.63482.43). [2] J. Robins. "A new approach to causal inference in mortality studies with a sustained exposure period-application to control of the healthy worker survivor effect". En. In: _"Math. Modelling"_ 7 (1986), pp. 1393-1512. [3] J. Kaufman. _Conceptos de Causalidad y Sesgos en Estudios Observacionales_. 2020. URL: [https://fivethirtyeight.com/features/science-isnt-broken/#part1](https://fivethirtyeight.com/features/science-isnt-broken/#part1) (visited on 07/09/2022). [4] A. Gelman and G. Imbens. _Why ask Why? Forward Causal Inference and Reverse Causal Questions_. Working Paper 19614. National Bureau of Economic Research, nov.. 2013. DOI: [10.3386/w19614](https://doi.org/10.3386%2Fw19614). URL: [http://www.nber.org/papers/w19614](http://www.nber.org/papers/w19614). [5] M. Ã. Hernán and J. M. Robins. _Causal Inference: What If_. Boca Raton: Chapman & Hall/CRC, 2020. URL: [https://cdn1.sph.harvard.edu/wp-content/uploads/sites/1268/2019/10/ci_hernanrobins_1oct19.pdf](https://cdn1.sph.harvard.edu/wp-content/uploads/sites/1268/2019/10/ci_hernanrobins_1oct19.pdf). [6] J. Pearl. "Statistics and Causal Inference: A Review". In: _TEST: An Official Journal of the Spanish Society of Statistics and Operations Research_ 12 (feb.. 2003), pp. 281-345. DOI: [10.1007/BF02595718](https://doi.org/10.1007%2FBF02595718). --- ## Referencias (2) . [7] A. Muriel, D. Hernández, and V. Abraira. "Modelos estructurales marginales: una herramienta útil que proporciona evidencia a los estudios observacionales". In: _Nefrologia_ 2.7 (2011), pp. 7-13. ISSN: 20137575. DOI: [10.3265/NefrologiaSuplementoExtraordinario.pre2011.Dec.11267](https://doi.org/10.3265%2FNefrologiaSuplementoExtraordinario.pre2011.Dec.11267). URL: [https://www.revistanefrologia.com/es-modelos-estructurales-marginales-una-herramienta-articulo-X2013757511000636](https://www.revistanefrologia.com/es-modelos-estructurales-marginales-una-herramienta-articulo-X2013757511000636). [8] C. Aschwanden. _Science Isnt Broken Its just a hell of a lot harder than we give it credit for_. 2015. URL: [https://fivethirtyeight.com/features/science-isnt-broken/#part1](https://fivethirtyeight.com/features/science-isnt-broken/#part1) (visited on 06/09/2022). [9] D. B. Rubin. "Randomization analysis of experimental data: The Fisher randomization test comment". In: _Journal of the American statistical association_ 75.371 (1980), pp. 591-593. [10] G. W. Imbens and D. B. Rubin. "Causality: The Basic Framework". In: _Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction_. Cambridge University Press, 2015, p. 3–22. DOI: [10.1017/CBO9781139025751.002](https://doi.org/10.1017%2FCBO9781139025751.002). [11] D. Tompsett, S. Vansteelandt, O. Dukes, et al. "gesttools: General Purpose G-Estimation in R". En. In: _Observational Studies_ 8.1 (2022), pp. 1-28. ISSN: 2767-3324. DOI: [10.1353/obs.2022.0003](https://doi.org/10.1353%2Fobs.2022.0003). URL: [https://muse.jhu.edu/article/856403](https://muse.jhu.edu/article/856403) (visited on 07/10/2022). [12] S. Vansteelandt and N. Keiding. "Invited commentary: G-computation-lost in translation?" En. In: _Am J Epidemiol_ 173.7 (mar.. 2011), pp. 739-742. --- ## Referencias (3) . [13] A. I. Naimi, S. R. Cole, and E. H. Kennedy. "An introduction to g methods". In: _International Journal of Epidemiology_ 46.2 (dic.. 2016), pp. 756-762. ISSN: 0300-5771. DOI: [10.1093/ije/dyw323](https://doi.org/10.1093%2Fije%2Fdyw323). eprint: https://academic.oup.com/ije/article-pdf/46/2/756/25421672/dyw323.pdf. URL: [https://doi.org/10.1093/ije/dyw323](https://doi.org/10.1093/ije/dyw323). [14] M. Bounthavong. _Using inverse probability of treatment weights & Marginal structural models to handle time-varying covariates_. 2018. URL: [https://rpubs.com/mbounthavong/IPTW_MSM_Tutorial](https://rpubs.com/mbounthavong/IPTW_MSM_Tutorial) (visited on 06/09/2022). [15] W. M. van der Wal and R. B. Geskus. "ipw: An R Package for Inverse Probability Weighting". In: _Journal of Statistical Software_ 43.13 (2011), pp. 1-23. URL: [http://www.jstatsoft.org/v43/i13/](http://www.jstatsoft.org/v43/i13/). [16] R. W. Platt, E. F. Schisterman, and S. R. Cole. "Time-modified Confounding". In: _American Journal of Epidemiology_ 170.6 (ago.. 2009), pp. 687-694. ISSN: 0002-9262. DOI: [10.1093/aje/kwp175](https://doi.org/10.1093%2Faje%2Fkwp175). eprint: https://academic.oup.com/aje/article-pdf/170/6/687/285709/kwp175.pdf. URL: [https://doi.org/10.1093/aje/kwp175](https://doi.org/10.1093/aje/kwp175).